package net.yegong.matrix;

import org.netlib.lapack.Sgeevx;
import org.netlib.lapack.Ssyevr;
import org.netlib.util.floatW;
import org.netlib.util.intW;

/**
 * @author cooper
 */
public class EigenvalueDecomposition {

	public EigenvalueDecomposition(SquareMatrix mat) {
		this(mat, false);
	}

	public EigenvalueDecomposition(SquareMatrix mat, boolean dOnly) {
		this.errors = MatrixHelper.eps();
		this.mat = mat;
		this.dOnly = dOnly;
	}

	public void setErrorsControl(float errors) {
		this.errors = errors;
	}

	/**
	 * The normal evd method. <div style="margin-left=4ex"> if
	 * <code>[V D] = evd(A)</code>, then there is <code>A * V = V * D</code>
	 * </div> You can get matrix V and matrix D through getV() or getD().
	 * Actually, you needn't invoke this method, getV() or getD() will do.
	 * 
	 * @see {@link getV()}, {@link getD()}, {@link getImageD()}
	 */
	public void calc() {
		dr = di = null;
		v = null;
		if (mat instanceof SymmetricMatrix) {
			SymmetricMatrix sy = (SymmetricMatrix) mat;
			calc(sy);
		}
		else {
			calc(mat);
		}
	}

	/**
	 * The partial evd method. You must use it like this. <div
	 * style="margin-left=4ex"><code>
	 * EigenvalueDecomposition evd = new EigenvalueDecomposition((SymmetricMatrix) m);<br />
	 * evd.calc(0, 1) // this compute the biggest eigenvalue & eigenvector<br />
	 * Matrix vectors = evd.getV();<br />
     * </code></div> You must remeber this, the non-symmetric can't take a
	 * partial evd. And the result is a normal evd.
	 * 
	 * @param indexFromLarge
	 * @param count
	 */
	public void partial(int indexFromLarge, int count) {
		dr = di = null;
		v = null;
		int upper = mat.nCol - indexFromLarge;
		int lower = upper - count + 1;
		if (upper > mat.nCol || lower < 0 || lower > upper)
			throw new IllegalArgumentException("The index of eigenvalues is illegal.");
		if (mat instanceof SymmetricMatrix) {
			SymmetricMatrix sy = (SymmetricMatrix) mat;
			calc(sy, lower, upper);
		}
		else {
			calc(mat);
		}
	}

	private void calc(SquareMatrix mat) {
		String jobvr = "N";
		int n = mat.nCol;
		float[] a = mat.getColumnMajorArrayCopy();
		float[] wr = new float[n];
		float[] wi = new float[n];
		float[] vl = null;
		float[] vr = null;
		intW ilo = new intW(0);
		intW ihi = new intW(0);
		float[] scale = new float[n];
		floatW abnrm = new floatW(0.0f);
		float[] rconde = new float[n];
		float[] rcondv = new float[n];
		float[] work = new float[1];
		int[] iwork = null;
		intW info = new intW(0);

		if (!dOnly) {
			jobvr = "V";
			vr = new float[n * n];
		}
		Sgeevx.sgeevx("B", "N", jobvr, "N", n, a, 0, n, wr, 0, wi, 0, vl, 0, n, vr, 0, n, ilo, ihi, scale, 0, abnrm, rconde, 0, rcondv, 0, work, 0, -1, iwork, 0, info);

		int lwork = (int) (work[0] + 0.1);
		work = new float[lwork];

		Sgeevx.sgeevx("B", "N", jobvr, "N", n, a, 0, n, wr, 0, wi, 0, vl, 0, n, vr, 0, n, ilo, ihi, scale, 0, abnrm, rconde, 0, rcondv, 0, work, 0, lwork, iwork, 0, info);

		if (!dOnly) {
			v = new SquareMatrix(vr, n, Matrix.FORTRAN);
		}
		dr = new DiagonalMatrix(wr);
		RowVector vec = new RowVector(wi);
		if (vec.maxAbsolute() >= vec.getElementsCount() * errors) {
			di = new DiagonalMatrix(wi);
		}
	}

	private void calc(SymmetricMatrix mat) {
		String jobz = "N";
		int n = mat.nCol;
		float[] a = mat.getColumnMajorArrayCopy();
		float abstol = errors;
		intW m = new intW(0);
		float[] w = new float[n];
		float[] z = null;
		int[] isuppz = null;
		float[] work = new float[1];
		int[] iwork = new int[1];
		intW info = new intW(0);

		if (!dOnly) {
			jobz = "V";
			z = new float[n * n];
			isuppz = new int[2 * n];
		}
		Ssyevr.ssyevr(jobz, "A", "U", n, a, 0, n, 0, 0, 0, 0, abstol, m, w, 0, z, 0, n, isuppz, 0, work, 0, -1, iwork, 0, -1, info);

		int lwork = (int) (work[0] + 0.1);
		int liwork = iwork[0];

		work = new float[lwork];
		iwork = new int[liwork];

		Ssyevr.ssyevr(jobz, "A", "U", n, a, 0, n, 0, 0, 0, 0, abstol, m, w, 0, z, 0, n, isuppz, 0, work, 0, lwork, iwork, 0, liwork, info);
		if (info.val < 0) {
			String s = "The " + -info.val + "-th parameter of lapack had an illegal value.";
			throw new RuntimeException(s);
		}
		if (!dOnly) {
			v = new SquareMatrix(z, n, Matrix.FORTRAN);
		}
		dr = new DiagonalMatrix(w);
		di = null;
	}

	private void calc(SymmetricMatrix mat, int lower, int upper) {
		int count = upper - lower + 1;

		String jobz = "N";
		int n = mat.nCol;
		float[] a = mat.getColumnMajorArrayCopy();
		float abstol = errors;
		intW m = new intW(0);
		float[] w = new float[n];
		float[] z = null;
		int[] isuppz = null;
		float[] work = new float[26 * n];
		int lwork = 26 * n;
		int[] iwork = new int[10 * n];
		int liwork = 10 * n;
		intW info = new intW(0);

		if (!dOnly) {
			jobz = "V";
			z = new float[n * count];
			isuppz = new int[2 * count];
		}
		Ssyevr.ssyevr(jobz, "I", "U", n, a, 0, n, 0, 0, lower, upper, abstol, m, w, 0, z, 0, n, isuppz, 0, work, 0, lwork, iwork, 0, liwork, info);
		if (info.val < 0) {
			String s = "The " + -info.val + "-th parameter of lapack had an illegal value.";
			throw new RuntimeException(s);
		}
		if (!dOnly) {
			v = new Matrix(z, n, count, Matrix.FORTRAN);
		}
		float[] w2 = new float[count];
		System.arraycopy(w, 0, w2, 0, count);
		dr = new DiagonalMatrix(w2);
		di = null;
	}

	public Matrix getV() {
		if (dr == null)
			calc();
		return v;
	}

	public DiagonalMatrix getD() {
		if (dr == null)
			calc();
		return dr;
	}

	public DiagonalMatrix getDImage() {
		if (mat instanceof SymmetricMatrix)
			return null;
		if (dr == null)
			calc();
		return di;
	}

	public boolean hasImage() {
		if (mat instanceof SymmetricMatrix)
			return false;
		if (dr == null)
			calc();
		return di != null;
	}

	protected final boolean dOnly;
	protected float errors;
	protected SquareMatrix mat;
	protected Matrix v;
	protected DiagonalMatrix dr;
	protected DiagonalMatrix di;
}
